clear all; close all; clc;
sigma = 13;
N = 1e6;
y_t = 0:0.1:65;
p_t = y_t/sigma^2 .* ...
exp(- y_t.^2 / 2 / sigma^2 );
x_n=rand(1, N);
y_n=sigma*sqrt(-2*log(x_n));
[h, y] = hist(y_n, 30);
inte = sum(h) * (y(2) - y(1));
p = h / inte;
figure(3);
plot(y_t, p_t, 'g', 'LineWidth', 5);
hold on; bar(y, p); hold off
xlabel('y'); ylabel('p(y)');
